##### This script runs the linear mixed-effects models reported in the Table OA4.

# The following set of models account for the possible BE and WE of country-round contextual variables
# However, the following set of models do not account for compositional effects

# Model A: two levels with RE at the COUNTRY-ROUND level
Model_A_CWC_ind_WE_BE <- lmer(red_sup ~ employed_CWC + unemployed_CWC + hinc_2_cop_CWC + hinc_3_dif_CWC + hinc_4_vdif_CWC + female_CWC + age_CWC + isced2_CWC + isced3_CWC + isced4_CWC + isced56_CWC + stfeco_CWC + lrscale_CWC + 
                        mean_LFPR_y0 + dev_mean_LFPR_y0 + mean_SOCEXP_y0 + dev_mean_SOCEXP_y0 + mean_GDP_y0 + dev_mean_GDP_y0 + mean_gini_disp_y0 + dev_mean_gini_disp_y0 + liberal + mediter + (1 | cntry_round), 
                      data = ess_4_ld_models, weights = pspwght)
summary(Model_A_CWC_ind_WE_BE)
check_collinearity(Model_A_CWC_ind_WE_BE)


# Model B: two levels with RE at the COUNTRY level
Model_B_CWC_ind_WE_BE <- lmer(red_sup ~ employed_CWC + unemployed_CWC + hinc_2_cop_CWC + hinc_3_dif_CWC + hinc_4_vdif_CWC + female_CWC + age_CWC + isced2_CWC + isced3_CWC + isced4_CWC + isced56_CWC + stfeco_CWC + lrscale_CWC + 
                        mean_LFPR_y0 + dev_mean_LFPR_y0 + mean_SOCEXP_y0 + dev_mean_SOCEXP_y0 + mean_GDP_y0 + dev_mean_GDP_y0 + mean_gini_disp_y0 + dev_mean_gini_disp_y0 + liberal + mediter + (1 | cntry), 
                      data = ess_4_ld_models, weights = pspwght)
summary(Model_B_CWC_ind_WE_BE)
check_collinearity(Model_B_CWC_ind_WE_BE)


# Model C: three levels with RE at ROUND level (3) and then COUNTRY-ROUND level (2)
Model_C_CWC_ind_WE_BE <- lmer(red_sup ~ employed_CWC + unemployed_CWC + hinc_2_cop_CWC + hinc_3_dif_CWC + hinc_4_vdif_CWC + female_CWC + age_CWC + isced2_CWC + isced3_CWC + isced4_CWC + isced56_CWC + stfeco_CWC + lrscale_CWC + 
                        mean_LFPR_y0 + dev_mean_LFPR_y0 + mean_SOCEXP_y0 + dev_mean_SOCEXP_y0 + mean_GDP_y0 + dev_mean_GDP_y0 + mean_gini_disp_y0 + dev_mean_gini_disp_y0 + liberal + mediter + (1 | cntry_round) + (1 | essround), 
                      data = ess_4_ld_models, weights = pspwght)
summary(Model_C_CWC_ind_WE_BE)
check_collinearity(Model_C_CWC_ind_WE_BE)


# Model D: three levels with RE at COUNTRY level (3) and then COUNTRY-ROUND level (2)
Model_D_CWC_ind_WE_BE <- lmer(red_sup ~ employed_CWC + unemployed_CWC + hinc_2_cop_CWC + hinc_3_dif_CWC + hinc_4_vdif_CWC + female_CWC + age_CWC + isced2_CWC + isced3_CWC + isced4_CWC + isced56_CWC + stfeco_CWC + lrscale_CWC + 
                        mean_LFPR_y0 + dev_mean_LFPR_y0 + mean_SOCEXP_y0 + dev_mean_SOCEXP_y0 + mean_GDP_y0 + dev_mean_GDP_y0 + mean_gini_disp_y0 + dev_mean_gini_disp_y0 + liberal + mediter + (1 | cntry_round) + (1 | cntry), 
                      data = ess_4_ld_models, weights = pspwght)
summary(Model_D_CWC_ind_WE_BE)
check_collinearity(Model_D_CWC_ind_WE_BE)


# Model E: cross-classified model
Model_E_CWC_ind_WE_BE <- lmer(red_sup ~ employed_CWC + unemployed_CWC + hinc_2_cop_CWC + hinc_3_dif_CWC + hinc_4_vdif_CWC + female_CWC + age_CWC + isced2_CWC + isced3_CWC + isced4_CWC + isced56_CWC + stfeco_CWC + lrscale_CWC + 
                        mean_LFPR_y0 + dev_mean_LFPR_y0 + mean_SOCEXP_y0 + dev_mean_SOCEXP_y0 + mean_GDP_y0 + dev_mean_GDP_y0 + mean_gini_disp_y0 + dev_mean_gini_disp_y0 + liberal + mediter + (1 | essround) + (1 | cntry), 
                      data = ess_4_ld_models, weights = pspwght)
summary(Model_E_CWC_ind_WE_BE)
check_collinearity(Model_E_CWC_ind_WE_BE)


# Model F: full model
Model_F_CWC_ind_WE_BE <- lmer(red_sup ~ employed_CWC + unemployed_CWC + hinc_2_cop_CWC + hinc_3_dif_CWC + hinc_4_vdif_CWC + female_CWC + age_CWC + isced2_CWC + isced3_CWC + isced4_CWC + isced56_CWC + stfeco_CWC + lrscale_CWC + 
                        mean_LFPR_y0 + dev_mean_LFPR_y0 + mean_SOCEXP_y0 + dev_mean_SOCEXP_y0 + mean_GDP_y0 + dev_mean_GDP_y0 + mean_gini_disp_y0 + dev_mean_gini_disp_y0 + liberal + mediter + (1 | cntry_round) + (1 | essround) + (1 | cntry), 
                      data = ess_4_ld_models, weights = pspwght)
summary(Model_F_CWC_ind_WE_BE)
check_collinearity(Model_F_CWC_ind_WE_BE)


# Model G: two levels with RE at the YEAR level
Model_G_CWC_ind_WE_BE <- lmer(red_sup ~ employed_CWC + unemployed_CWC + hinc_2_cop_CWC + hinc_3_dif_CWC + hinc_4_vdif_CWC + female_CWC + age_CWC + isced2_CWC + isced3_CWC + isced4_CWC + isced56_CWC + stfeco_CWC + lrscale_CWC + 
                        mean_LFPR_y0 + dev_mean_LFPR_y0 + mean_SOCEXP_y0 + dev_mean_SOCEXP_y0 + mean_GDP_y0 + dev_mean_GDP_y0 + mean_gini_disp_y0 + dev_mean_gini_disp_y0 + liberal + mediter + (1 | essround), 
                      data = ess_4_ld_models, weights = pspwght)
summary(Model_G_CWC_ind_WE_BE)
check_collinearity(Model_G_CWC_ind_WE_BE)


# Plotting the model results
tab_model(Model_A_CWC_ind_WE_BE, Model_B_CWC_ind_WE_BE, Model_C_CWC_ind_WE_BE, Model_D_CWC_ind_WE_BE, Model_E_CWC_ind_WE_BE, Model_F_CWC_ind_WE_BE, Model_G_CWC_ind_WE_BE, digits = 3, digits.p = 4, digits.re = 4)



#################################################################################
# Storing and exporting results of the seven models into a Table OA4 of the paper #
#################################################################################

# Creating an empty data.frame into which the results will be stored
Models_CWC_ind_WE_BE <- data.frame(c("Intercept", "", "employed_CWC", "", "unemployed_CWC", "", 
                             "hinc_2_cop_CWC", "", "hinc_3_dif_CWC", "", "hinc_4_vdif_CWC", "", 
                             "female_CWC", "", "age_CWC", "", "isced2_CWC", "", "isced3_CWC", "", 
                             "isced4_CWC", "", "isced56_CWC", "", "stfeco_CWC", "", "lrscale_CWC", "", 
                             "mean_LFPR_y0", "", "dev_mean_LFPR_y0", "", "mean_SOCEXP_y0", "", "dev_mean_SOCEXP_y0", "", 
                             "mean_GDP_y0", "", "dev_mean_GDP_y0", "", "mean_gini_disp_y0", "", "dev_mean_gini_disp_y0", "", 
                             "liberal", "", "mediter", "", "Var_cntry", "Var_essround", 
                             "Var_cntry_round", "Var_individual", "ICC_cntry", "ICC_essround", "ICC_cntry_round"), matrix(data = rep(NA, 55*7), nrow = 55))
colnames(Models_CWC_ind_WE_BE) <- c("Parameter", "Model A", "Model B", "Model C", "Model D", "Model E", "Model F", "Model G")

# Using the programmed function to create Table OA4 of the article
Models_CWC_ind_WE_BE <- assign_fe_results(model_name = Model_A_CWC_ind_WE_BE, assignment_object = Models_CWC_ind_WE_BE, column_number = 2)
Models_CWC_ind_WE_BE <- assign_fe_results(model_name = Model_B_CWC_ind_WE_BE, assignment_object = Models_CWC_ind_WE_BE, column_number = 3)
Models_CWC_ind_WE_BE <- assign_fe_results(model_name = Model_C_CWC_ind_WE_BE, assignment_object = Models_CWC_ind_WE_BE, column_number = 4)
Models_CWC_ind_WE_BE <- assign_fe_results(model_name = Model_D_CWC_ind_WE_BE, assignment_object = Models_CWC_ind_WE_BE, column_number = 5)
Models_CWC_ind_WE_BE <- assign_fe_results(model_name = Model_E_CWC_ind_WE_BE, assignment_object = Models_CWC_ind_WE_BE, column_number = 6)
Models_CWC_ind_WE_BE <- assign_fe_results(model_name = Model_F_CWC_ind_WE_BE, assignment_object = Models_CWC_ind_WE_BE, column_number = 7)
Models_CWC_ind_WE_BE <- assign_fe_results(model_name = Model_G_CWC_ind_WE_BE, assignment_object = Models_CWC_ind_WE_BE, column_number = 8)
View(Models_CWC_ind_WE_BE)

### Manual assignment of variance components and ICCs
# Model A: RE at cntry_round
Models_CWC_ind_WE_BE[Models_CWC_ind_WE_BE$Parameter == "Var_cntry_round", "Model A"] <- as.data.frame(VarCorr(Model_A_CWC_ind_WE_BE))[, "vcov"][1]
Models_CWC_ind_WE_BE[Models_CWC_ind_WE_BE$Parameter == "Var_individual", "Model A"] <- as.data.frame(VarCorr(Model_A_CWC_ind_WE_BE))[, "vcov"][2]
icc(Model_A_CWC_ind_WE_BE, by_group = TRUE)
Models_CWC_ind_WE_BE[Models_CWC_ind_WE_BE$Parameter == "ICC_cntry_round", "Model A"] <- (as.data.frame(VarCorr(Model_A_CWC_ind_WE_BE))[, "vcov"][1]/sum(as.data.frame(VarCorr(Model_A_CWC_ind_WE_BE))[, "vcov"]))*100

# Model B: RE at cntry
Models_CWC_ind_WE_BE[Models_CWC_ind_WE_BE$Parameter == "Var_cntry", "Model B"] <- as.data.frame(VarCorr(Model_B_CWC_ind_WE_BE))[, "vcov"][1]
Models_CWC_ind_WE_BE[Models_CWC_ind_WE_BE$Parameter == "Var_individual", "Model B"] <- as.data.frame(VarCorr(Model_B_CWC_ind_WE_BE))[, "vcov"][2]
icc(Model_B_CWC_ind_WE_BE, by_group = TRUE)
Models_CWC_ind_WE_BE[Models_CWC_ind_WE_BE$Parameter == "ICC_cntry", "Model B"] <- (as.data.frame(VarCorr(Model_B_CWC_ind_WE_BE))[, "vcov"][1]/sum(as.data.frame(VarCorr(Model_B_CWC_ind_WE_BE))[, "vcov"]))*100

# Model C: RE at essround and cntry_round
Models_CWC_ind_WE_BE[Models_CWC_ind_WE_BE$Parameter == "Var_essround", "Model C"] <- as.data.frame(VarCorr(Model_C_CWC_ind_WE_BE))[, "vcov"][2]
Models_CWC_ind_WE_BE[Models_CWC_ind_WE_BE$Parameter == "Var_cntry_round", "Model C"] <- as.data.frame(VarCorr(Model_C_CWC_ind_WE_BE))[, "vcov"][1]
Models_CWC_ind_WE_BE[Models_CWC_ind_WE_BE$Parameter == "Var_individual", "Model C"] <- as.data.frame(VarCorr(Model_C_CWC_ind_WE_BE))[, "vcov"][3]
icc(Model_C_CWC_ind_WE_BE, by_group = TRUE)
Models_CWC_ind_WE_BE[Models_CWC_ind_WE_BE$Parameter == "ICC_essround", "Model C"] <- (as.data.frame(VarCorr(Model_C_CWC_ind_WE_BE))[, "vcov"][2]/sum(as.data.frame(VarCorr(Model_C_CWC_ind_WE_BE))[, "vcov"]))*100
Models_CWC_ind_WE_BE[Models_CWC_ind_WE_BE$Parameter == "ICC_cntry_round", "Model C"] <- (as.data.frame(VarCorr(Model_C_CWC_ind_WE_BE))[, "vcov"][1]/sum(as.data.frame(VarCorr(Model_C_CWC_ind_WE_BE))[, "vcov"]))*100

# Model D: three levels with RE at COUNTRY level (3) and then COUNTRY-ROUND level (2)
Models_CWC_ind_WE_BE[Models_CWC_ind_WE_BE$Parameter == "Var_cntry", "Model D"] <- as.data.frame(VarCorr(Model_D_CWC_ind_WE_BE))[, "vcov"][2]
Models_CWC_ind_WE_BE[Models_CWC_ind_WE_BE$Parameter == "Var_cntry_round", "Model D"] <- as.data.frame(VarCorr(Model_D_CWC_ind_WE_BE))[, "vcov"][1]
Models_CWC_ind_WE_BE[Models_CWC_ind_WE_BE$Parameter == "Var_individual", "Model D"] <- as.data.frame(VarCorr(Model_D_CWC_ind_WE_BE))[, "vcov"][3]
icc(Model_D_CWC_ind_WE_BE, by_group = TRUE)
Models_CWC_ind_WE_BE[Models_CWC_ind_WE_BE$Parameter == "ICC_cntry", "Model D"] <- (as.data.frame(VarCorr(Model_D_CWC_ind_WE_BE))[, "vcov"][2]/sum(as.data.frame(VarCorr(Model_D_CWC_ind_WE_BE))[, "vcov"]))*100
Models_CWC_ind_WE_BE[Models_CWC_ind_WE_BE$Parameter == "ICC_cntry_round", "Model D"] <- (as.data.frame(VarCorr(Model_D_CWC_ind_WE_BE))[, "vcov"][1]/sum(as.data.frame(VarCorr(Model_D_CWC_ind_WE_BE))[, "vcov"]))*100

# Model E: RE at country and ess_round
Models_CWC_ind_WE_BE[Models_CWC_ind_WE_BE$Parameter == "Var_cntry", "Model E"] <- as.data.frame(VarCorr(Model_E_CWC_ind_WE_BE))[, "vcov"][1]
Models_CWC_ind_WE_BE[Models_CWC_ind_WE_BE$Parameter == "Var_essround", "Model E"] <- as.data.frame(VarCorr(Model_E_CWC_ind_WE_BE))[, "vcov"][2]
Models_CWC_ind_WE_BE[Models_CWC_ind_WE_BE$Parameter == "Var_individual", "Model E"] <- as.data.frame(VarCorr(Model_E_CWC_ind_WE_BE))[, "vcov"][3]
icc(Model_E_CWC_ind_WE_BE, by_group = TRUE)
Models_CWC_ind_WE_BE[Models_CWC_ind_WE_BE$Parameter == "ICC_cntry", "Model E"] <- (as.data.frame(VarCorr(Model_E_CWC_ind_WE_BE))[, "vcov"][1]/sum(as.data.frame(VarCorr(Model_E_CWC_ind_WE_BE))[, "vcov"]))*100
Models_CWC_ind_WE_BE[Models_CWC_ind_WE_BE$Parameter == "ICC_essround", "Model E"] <- (as.data.frame(VarCorr(Model_E_CWC_ind_WE_BE))[, "vcov"][2]/sum(as.data.frame(VarCorr(Model_E_CWC_ind_WE_BE))[, "vcov"]))*100

# Model F: RE at all three effects
Models_CWC_ind_WE_BE[Models_CWC_ind_WE_BE$Parameter == "Var_cntry", "Model F"] <- as.data.frame(VarCorr(Model_F_CWC_ind_WE_BE))[, "vcov"][2]
Models_CWC_ind_WE_BE[Models_CWC_ind_WE_BE$Parameter == "Var_essround", "Model F"] <- as.data.frame(VarCorr(Model_F_CWC_ind_WE_BE))[, "vcov"][3]
Models_CWC_ind_WE_BE[Models_CWC_ind_WE_BE$Parameter == "Var_cntry_round", "Model F"] <- as.data.frame(VarCorr(Model_F_CWC_ind_WE_BE))[, "vcov"][1]
Models_CWC_ind_WE_BE[Models_CWC_ind_WE_BE$Parameter == "Var_individual", "Model F"] <- as.data.frame(VarCorr(Model_F_CWC_ind_WE_BE))[, "vcov"][4]
icc(Model_F_CWC_ind_WE_BE, by_group = TRUE)
Models_CWC_ind_WE_BE[Models_CWC_ind_WE_BE$Parameter == "ICC_cntry", "Model F"] <- (as.data.frame(VarCorr(Model_F_CWC_ind_WE_BE))[, "vcov"][2]/sum(as.data.frame(VarCorr(Model_F_CWC_ind_WE_BE))[, "vcov"]))*100
Models_CWC_ind_WE_BE[Models_CWC_ind_WE_BE$Parameter == "ICC_essround", "Model F"] <- (as.data.frame(VarCorr(Model_F_CWC_ind_WE_BE))[, "vcov"][3]/sum(as.data.frame(VarCorr(Model_F_CWC_ind_WE_BE))[, "vcov"]))*100
Models_CWC_ind_WE_BE[Models_CWC_ind_WE_BE$Parameter == "ICC_cntry_round", "Model F"] <- (as.data.frame(VarCorr(Model_F_CWC_ind_WE_BE))[, "vcov"][1]/sum(as.data.frame(VarCorr(Model_F_CWC_ind_WE_BE))[, "vcov"]))*100

# Model G: RE at essround
Models_CWC_ind_WE_BE[Models_CWC_ind_WE_BE$Parameter == "Var_essround", "Model G"] <- as.data.frame(VarCorr(Model_G_CWC_ind_WE_BE))[, "vcov"][1]
Models_CWC_ind_WE_BE[Models_CWC_ind_WE_BE$Parameter == "Var_individual", "Model G"] <- as.data.frame(VarCorr(Model_G_CWC_ind_WE_BE))[, "vcov"][2]
icc(Model_G_CWC_ind_WE_BE, by_group = TRUE)
Models_CWC_ind_WE_BE[Models_CWC_ind_WE_BE$Parameter == "ICC_essround", "Model G"] <- (as.data.frame(VarCorr(Model_G_CWC_ind_WE_BE))[, "vcov"][1]/sum(as.data.frame(VarCorr(Model_G_CWC_ind_WE_BE))[, "vcov"]))*100

# Exporting Table OA4 to a csv file
write.csv(x = Models_CWC_ind_WE_BE, file = "table_OA4_Models_CWC_ind_BEs_WEs.csv")



### Collecting the coefficient estimates for the LFPR BE contextual variable (when individual-level variables are CWC) and storing them to a data.frame
LFPR_CWC_ind_BE <- data.frame(c("Model A", "Model B", "Model C", "Model D", "Model E", "Model F", "Model G"),
                      matrix(data = rep(NA, 7*3), nrow = 7))
colnames(LFPR_CWC_ind_BE) <- c("Model type", "Estimate", "Std. Error", "df")

LFPR_CWC_ind_BE[LFPR_CWC_ind_BE$`Model type` == "Model A", 2:4] <- coef(summary(Model_A_CWC_ind_WE_BE))[15, c("Estimate", "Std. Error", "df")]
LFPR_CWC_ind_BE[LFPR_CWC_ind_BE$`Model type` == "Model B", 2:4] <- coef(summary(Model_B_CWC_ind_WE_BE))[15, c("Estimate", "Std. Error", "df")]
LFPR_CWC_ind_BE[LFPR_CWC_ind_BE$`Model type` == "Model C", 2:4] <- coef(summary(Model_C_CWC_ind_WE_BE))[15, c("Estimate", "Std. Error", "df")]
LFPR_CWC_ind_BE[LFPR_CWC_ind_BE$`Model type` == "Model D", 2:4] <- coef(summary(Model_D_CWC_ind_WE_BE))[15, c("Estimate", "Std. Error", "df")]
LFPR_CWC_ind_BE[LFPR_CWC_ind_BE$`Model type` == "Model E", 2:4] <- coef(summary(Model_E_CWC_ind_WE_BE))[15, c("Estimate", "Std. Error", "df")]
LFPR_CWC_ind_BE[LFPR_CWC_ind_BE$`Model type` == "Model F", 2:4] <- coef(summary(Model_F_CWC_ind_WE_BE))[15, c("Estimate", "Std. Error", "df")]
LFPR_CWC_ind_BE[LFPR_CWC_ind_BE$`Model type` == "Model G", 2:4] <- coef(summary(Model_G_CWC_ind_WE_BE))[15, c("Estimate", "Std. Error", "df")]

# Computing the 95% confidence interval (based on student?s t-distribution)
LFPR_CWC_ind_BE$lower_bound <- LFPR_CWC_ind_BE$Estimate - LFPR_CWC_ind_BE$`Std. Error`*qt(p = 0.975, df = LFPR_CWC_ind_BE$df)
LFPR_CWC_ind_BE$upper_bound <- LFPR_CWC_ind_BE$Estimate + LFPR_CWC_ind_BE$`Std. Error`*qt(p = 0.975, df = LFPR_CWC_ind_BE$df)


### Collecting the coefficient estimates for the LFPR WE contextual variable (when individual-level variables are CWC) and storing them to a data.frame
LFPR_CWC_ind_WE <- data.frame(c("Model A", "Model B", "Model C", "Model D", "Model E", "Model F", "Model G"),
                      matrix(data = rep(NA, 7*3), nrow = 7))
colnames(LFPR_CWC_ind_WE) <- c("Model type", "Estimate", "Std. Error", "df")

LFPR_CWC_ind_WE[LFPR_CWC_ind_WE$`Model type` == "Model A", 2:4] <- coef(summary(Model_A_CWC_ind_WE_BE))[16, c("Estimate", "Std. Error", "df")]
LFPR_CWC_ind_WE[LFPR_CWC_ind_WE$`Model type` == "Model B", 2:4] <- coef(summary(Model_B_CWC_ind_WE_BE))[16, c("Estimate", "Std. Error", "df")]
LFPR_CWC_ind_WE[LFPR_CWC_ind_WE$`Model type` == "Model C", 2:4] <- coef(summary(Model_C_CWC_ind_WE_BE))[16, c("Estimate", "Std. Error", "df")]
LFPR_CWC_ind_WE[LFPR_CWC_ind_WE$`Model type` == "Model D", 2:4] <- coef(summary(Model_D_CWC_ind_WE_BE))[16, c("Estimate", "Std. Error", "df")]
LFPR_CWC_ind_WE[LFPR_CWC_ind_WE$`Model type` == "Model E", 2:4] <- coef(summary(Model_E_CWC_ind_WE_BE))[16, c("Estimate", "Std. Error", "df")]
LFPR_CWC_ind_WE[LFPR_CWC_ind_WE$`Model type` == "Model F", 2:4] <- coef(summary(Model_F_CWC_ind_WE_BE))[16, c("Estimate", "Std. Error", "df")]
LFPR_CWC_ind_WE[LFPR_CWC_ind_WE$`Model type` == "Model G", 2:4] <- coef(summary(Model_G_CWC_ind_WE_BE))[16, c("Estimate", "Std. Error", "df")]

# Computing the 95% confidence interval (based on student?s t-distribution)
LFPR_CWC_ind_WE$lower_bound <- LFPR_CWC_ind_WE$Estimate - LFPR_CWC_ind_WE$`Std. Error`*qt(p = 0.975, df = LFPR_CWC_ind_WE$df)
LFPR_CWC_ind_WE$upper_bound <- LFPR_CWC_ind_WE$Estimate + LFPR_CWC_ind_WE$`Std. Error`*qt(p = 0.975, df = LFPR_CWC_ind_WE$df)



# Deleting model results from the workspace (as they are no longer needed)
rm(Model_A_CWC_ind_WE_BE, Model_B_CWC_ind_WE_BE, Model_C_CWC_ind_WE_BE, Model_D_CWC_ind_WE_BE, Model_E_CWC_ind_WE_BE, Model_F_CWC_ind_WE_BE, Model_G_CWC_ind_WE_BE)
